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Abstract 

We investigate dynamics of an inference algorithm termed the belief propagation (BP) when 
employed in spin glass (SG) models and show that its macroscopic behaviors can be traced by 
recursive updates of certain auxiliary field distributions whose stationary state reproduces the 
replica symmetric solution offered by the equilibrium analysis. We further provide a compact 
expression for the instability condition of the BP's fixed point which turns out to be identical to 
that of instability for breaking the replica symmetry in equilibrium when the number of couplings 
per spin is infinite. This correspondence is extended to a SG model of finite connectivity to 
determine the phase diagram, which is numerically supported. 

PACS numbers: 89.90.+n, 02.50.-r, 05.50.+q, 75.10.Hk 
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Recently, there is growing interest in a similarity between researches on spin glass (SG) 
and information processing (IP) Since employment of methods from SG theory provided 
significant progresses for several problems related to IP such as machine learning 0, error- 
correcting [[| [|, |5| and spreading codes |7|, it is natural to expect that the opposite 
direction might be possible. 

The purpose of this article is to show such an example. More specifically, we show that 
investigating dynamics of an iterative inference algorithm termed the belief propagation (BP) 
which has been developed in IP research [^, 10 1 when employed in SG models provides a new 



understanding about thermodynamical properties of SG. We show that the replica symmetric 
(RS) solution known in the equilibrium analysis can be characterized as a macro scopically 
stationary state in BP. We also provide a compact expression of the microscopic instability 
condition around the fixed point in the BP dynamics which turns out to be identical to that 
of instability for breaking the replica symmetry in equilibrium termed the Almeida-Thouless 



(AT) instability [ I 1 1 when the number of connectivity per spin is infinite. Efficacy of this 
expression for a sparsely connected SG model is also numerically supported. 
We here take up a family of Ising SG models define by Hamiltonian 

M 

h(s\j) = -j:j» n a) 

where £(/i) denotes a set of indices which are connected to a quenched coupling J M . We 
assume that each coupling is independently generated from an identical distribution 



We further assume that for each /i, £(/i) is composed of randomly selected K ~ 0(1) spin 
indices and each spin index I is concerned with C couplings the set of which is denoted as 
Ai(l). Jo > and J > are parameters to control the mean and the standard deviation of 
J M , respectively, which naturally links the current system (P to the Sherrington-Kirkpatrick 
(SK) model [Q in the case of K = 2 and C ~ O(N) and to sparsely connected SG models 



|T3| , p^| , |15| in general. 

A major goal of statistical mechanics in the current system is to calculate the mi- 
croscopic spin average mi = Tig Si exp [— (3H(S\ J)] /Trgexp [— f3H(S\ J)} from given 
Hamiltonian ([I]). This is formally identical to an inference problem for a posterior dis- 
tribution P(S\J) oc Il^Li P(J/j,\S) derived from a conditional probability P(J^\S) = 



exp \J3 J M Uiec(n) Si\ /Ej M=± j /v ^exp [PJ^Ui&c^) S i\ and a uniform prior, which can be ex- 
pressed in a bipartite graph as Figure [I] (a). In this expression, spins and couplings are 
denoted as two different types of nodes and are linked by edges when they are directly 
connected, which is useful to explicitly represent statistical dependences between estimation 
variables (spins) and observed data (couplings). 

BP is an iterative algorithm defined over the bipartite graph to calculate the spin average 
for a given set of couplings J = (J M ) || [TD[. In the current system, this is performed by 



passing beliefs (or messages) between the two types of nodes via edges at each update as 



'+ 1 ~- tanh/3J M J] m U> ( 3 ) 

keC((i)\l 



mjj = tanh [ J2 tanh 1 ™li ] > ( 4 ) 



vfeM(Q\/x 

where beliefs m*i and m*j are parameters to represent auxiliary distributions at tth update 
as P(S t \{J^}) = (l+m*,5i)/2 and P{J,\S h {J u ^}) = Tr 8m P(J»\S)P(S\{J^}) oc 
(1 + m* l Si)/2, respectively. stands for a set of spin indices which belong to 

other than / and similarly to Ai(l)\fi. Calculating iteratively, the estimate of the spin 
average at tth update is provided as 



VTL, 



tanh ^2 tanh 1 m* d . (5) 



It is known that BP provides the exact spin average by the convergent solution when the 
bipartite graph is free from cycles (Figure [l] (b)). Actually, BP is a very similar scheme to the 
transfer matrix method (TMM) or the Bethe approximation [lq, 17] which is frequently used 



in physics and the current statement can be regarded as a generalization of a known property 
of TMM that offers the exact results for a one dimensional lattice or a tree. However, BP 
still has a possibility to introduce something new into physics since it is explicitly expressed 
as an algorithm and such view point has been rare in the research on matters. This strongly 
motivates us to examine its dynamical properties, which we will focus on hereafter. 

Let us first discuss the macroscopic behavior of the BP dynamics (|3p and (||). Although 
the current randomly constructed system is not free from cycles, it can be shown that the 
typical length of the cycles grows as 0(ln N) with respect to the system size N as long as C is 
0(1) [|T^], which implies that the self-interaction from the past state is presumably negligible 
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in the thermodynamic limit. On the other hand, the self- interaction is also expected as suffi- 
ciently small even if C is large since the strength of the coupling becomes weak as 0(C~ 1//2 ). 
This and eqs. (|3|) and (£|) imply that the time evolution of the macroscopic distributions of 
beliefs tt'(x) = (1/NC) Z? =1 E^M(i) S(x-m%) and = (1/NC) E^i E^( ^x-m^) 
is likely to be well captured by recursive equations 



7T 



K-l I / K-l 

t+1 



x) = J _Q dxrfixi) [6 \x - tanh/3J _Q x,j y , (6) 
7r*(x) = J Yl dx^^x^d I x — tanh I ^ tanh -1 x^ J J , (7) 



where (• • ■) j represents the average with respect to J following distribution @. 

The validity of the current argument and its link to the replica symmetric (RS) ansatz 
in the equilibrium analysis have been shown already for finite C 0, ^0|. Here, we further 



show that these can be extended to the case of infinite C even if the AT stability of the RS 
solution is broken in equilibrium. 

When C becomes infinite, it is more convenient to deal with an auxiliary field of finite 
strength = T,»eM(i)\n tanh -1 m*, ~ Ei/eM(Q\f* ™li rather than m*, since becomes 
infinitesimal. Due to the central limit theorem, the distribution of the auxiliary field = 
(1/NC) Y4L1 Z)|teM(0 ~~ fyftf) can ^ e re S ar ded as a Gaussian 

y/2%F 



P*W = / II dx^ix^h - tanh" 1 ^) ~ — ==exp - v '" OF/ 7 , (8) 



where E 1 and are the average and the variance to parameterize the Gaussian distribution 
p t (h), respectively. The expression in the middle implies 7r*(x) = / dhp t (h)5(x — tanh(/i)). 
Plugging this into eq. @ and recursively employing eq. we obtain a compact expression 
for the update of E t and F l as 

E t+1 = PJ (M 1 )*' 1 , F t+1 = (3 2 J 2 (Q'Y' 1 , (9) 
M* = J Dztanh(VF~ t z + E 1 ), Q l = J Dz tanh 2 (VF*z + E*), (10) 

where = exp[— z 2 /2]/y / 27r and M* and can be expressed as M* ~ (1/N) J2iLi — 
(1/N) Eti m \ and Q* ~ (1/iV) E^m^) 2 ^ (1/N) E^i(W) 2 , respectively, due to the law 
of large numbers. Eqs. (§) and (|10|) serve as alternatives of eqs. (H) and (|7]). 

It should be noticed here that these equations can be regarded as the forward iteration 
of the saddle point equations to obtain the RS solution in the replica analysis of the multi- 
spin interaction infinite connectivity SG models ffl and, in particular, of the SK model for 



K = 2 [12]. In order to confirm the validity of the above argument, we compared the time 
evolution of the belief update (g) and (|) (BP) with that of eqs. (g) and (0) (RS) for the 
SK (K = 2) model, which is shown in Figures |2] (a) and (b). We also compared them with 
the trajectory of the naive iteration of the BP's fixed point condition 

m,=tanh| V (3J^ J[ m k - W E i II m A (1 - ™?H I > (H) 



f e ^ n e (w 2 e ( n J 

y/iex(0 keC(n)\l h<eM(1) jec{fi)\i Ve£(M)VJ / 



(TAP) which can be obtained inserting m^j ~ mi — (1 — m^)m^i to the fixed point of eqs. 
(§) and (H) m*j = m^, = m M / and = m;. This becomes identical to the famous 



Thouless-Anderson-Palmer (TAP) equation of the SK model, in particular, for K = 2 pT 



The experiments were performed for J = 1.5, 0.5 keeping J = 1 and T = 0.5, where the 
AT stability of the RS solution in equilibrium is satisfied for J — 1-5 but broken for J — 0.5 
|TT| . Figures [| (a) and (b) show that BP and RS exhibit excellent consistency with respect 
to the macroscopic variables irrespectively of whether the AT stability is satisfied or not. 
This strongly validates the reduction from BP (|3]) and (||) to the macroscopic dynamics (0) 
and (|TUD. On the other hand, TAP is considerably different from the others. In a sense, 
this may be natural because naively iterating eq. (|TT]) is just one of procedures for obtaining 
a solution and its trajectory in dynamics does not necessarily have any consistency with 
BP or RS while the BP's fixed point is correctly characterized by the TAP equation (|TTD 



which does have a certain relation to the RS solution in equilibrium as shown in [22] . These 
figures also imply that the dynamics of BP cannot be traced by a closed set of equations 
with respect to singly indexed variables m\ even for C — > oo while the fixed point condition 
in this limit is provided as coupled equations of mi (|TT|), which is also observed in a similar 
system [Q. 

Although Figures |2| (a) and (b) show that the macroscopic variables rapidly converge to 
those of the RS solution in BP, this does not imply that BP microscopically converges to a 
certain solution. In order to probe this microscopic convergence, we examined the squared 
difference of spin averages between successive updates D l = (l/N) Y4=i{ m \ ~ m « _1 ) 2 ) ^ ne 
time evolution of which is shown in the insets of Figures |2] (a) and (b). These illustrate 
that the (microscopic) local stability of the BP's fixed point can be broken even if the 
macroscopic behavior seems to converge, which cannot be detected by only examining the 
reduced macroscopic dynamics (§) and (jlOD . 
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In order to characterize such instability, we next turn to the stability analysis of the 
BP updates and Linearizing the updates with respect to the auxiliary field h^i = 
tanh -1 m^x around a fixed point solution = m^, we obtain a dynamics of the auxiliary 



field fluctuation Sh^ as 



tanh f3J u l\kec(u)\i^uk ^ 1 



m 



2 



8hp= £ —r^^^^ E ( 12 ) 

ueM(i)\v l - (t&nh (3, J v Y\kec(v)\i m vk) jec(u)\i m ^ 



Analytically solving this linearized equation for a large graph is generally difficult. However, 
since the current system is randomly constructed, the self-interaction of from the past 
can be considered as small as those of beliefs are. This implies that the time evolution of 
the fluctuation distribution f l (y) = (1/NC) J2iLi H^m{i) ~~ ^Kii) can t> e provided by a 
functional equation 

f t+ \y) = n T[dy,if t {y,i) 

n=i 1=1 

I A V tanh^nf^ 1 ^ y^ j-^j ^\ 

x \ 6 \y - £ — / — K r \ 2 x £ ~~ x ^ ) > ( 13 ) 

\ ^ m=i 1- (tanh^nfji 1 ^) 1=1 ^ /'.y 

where (■ • •) ^ a ; denotes the average over and x^i following eq. (0) and the stationary 
distribution of tt*(x) = vr(x), respectively, and the stability of the BP's fixed point can be 
characterized by whether the stationary solution f t (y) = f(y) = 8{y) is stable or not in 
update (|I3D. This formulation makes analytical investigation possible to a certain extent. 

In order to connect eq. ([13]) to the existing analysis, let us first investigate the limit 
C — > oo for which much more results are known compared to the case of finite C. 
Due to the central limit theorem, the distribution of the field fluctuation can be as- 



sumed as a Gaussian f t {y) = (1/v27tP) exp [—(y — a 1 ) 2 / (26*)], where a* and h l are the 
mean and the variance of the distribution, respectively. Plugging this expression into 
eq. ([Hp offers update rules with respect to a* and 6* as a t+1 = /3J M K ~ 2 (1 — Q)a f and 
b t+l = (p,J) 2 Q K - 2 JDz(l- tanh 2 (V^ + E)f (6* + (a*) 2 ), where M, Q, E and F repre- 
sent the convergent solutions of eqs. (§) and fllQ|) . In order to examine the stability of 
f(y) = 5(y), we linearize these equations around a* = b l = 0, which provides the critical 
condition of the instability with respect to the growth of h l 

{/3J) 2 Q K - 2 J Dz (l - tanh 2 (^Fz + E)f = 1, (14) 
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which becomes identical to that of the AT stability for the infinite range multi-spin inter- 



action SG models and, in particular, for the SK model when K = 2 flllfl . Furthermore, 
in the case of the SK model (K = 2), the critical condition with respect to a 1 around the 
paramagnetic solution M = Q = corresponds to the para-ferromagnetic transition. These 
mean that the two different phase transitions from the paramagnetic solution can be linked 
in a unified framework to the dynamical instabilities of BP by eq. (|T3"D. 

When C is finite, one can numerically perform the stability analysis employing eq. (|L3|) . 
the detail of which will be reported elsewhere. In addition, analytical investigation becomes 
possible for K = 2 as follows since transitions from the paramagnetic solution in this case 
occur due to the local instability. 

For a small /3, the paramagnetic solution ir(x) = n(x) = S(x) {m^ = rh^i = 
0) expresses the correct stable fixed point of the BP dynamics. Inserting this into 
eq. (|13|) does not provide a closed set of equations with respect to a finite num- 
ber of parameters since f t {y) is no more a Gaussian. However, assuming — 
S(y), the stability analysis can be reduced to coupled equations with respect to the 
mean and the variance of f t {y) as a t+1 = [C — 1) (tanhf3J')ja t and b t+1 = (C — 



1) (^tanh 2 (3 J) j 6* + (^tanh 2 f3J) j - (ta,nh(3J}j) (a*) 2 J. Linearizing these around a 1 = 
b l = provides the critical conditions with respect to the growth of a 1 and b % as 

(C -l)(tanh/3J)j = 1, (15) 
(C - 1) (tanh 2 j = 1, (16) 

respectively. It should be mentioned that a similar condition to eq. ([16]) was once obtained 
for a SG model on the Bethe lattice |53[ while eq. (|1|) was not. However, the current 
scheme may be superior to that employed in 23] as the expression (p~3| ) is compact and, 



therefore, can be easily extended to the case of multi-spin interaction (K > 3) with the aid 
of numerical methods while such extension requires higher order perturbation and becomes 
highly complicated in the other scheme. 

Eqs. fll5|) and (]16|) might correspond to the para-ferromagnetic and the para-SG phase 
transitions, respectively, since they do in the limit C — ► oo. In order to examine this, we 
performed numerical experiments for N = 2000 and C — 4. Although further investigation 
may be necessary to prove correctness, the data obtained from 100 experiments of 20000 
Monte Carlo steps per spin exhibit good consistency with the analytical expressions ([T5|) 



and ([TBI) indicating that the correspondence between the phase transitions in equilibrium 
and the dynamical instabilities of BP holds for finite C as well (Figure |3|). 

In summary, we have investigated dynamical behavior of BP when employed in SG mod- 
els. We have shown that the time evolution of macroscopic variables can be well captured 
by recursive updates of auxiliary field distributions which becomes identical to the forward 
iteration of the saddle point equations under the RS ansatz in the replica analysis. We have 
further shown that the dynamical instability of the BP's fixed point is closely related to the 
AT instability of the RS solution, which has been numerically supported. 

Relationship between the current scheme and an existing AT analysis for finite connec- 



tivity SG models [24] that generally requires complicated calculation and is not frequently 
employed in practice is under investigation. Besides this, extension of the current framework 



to the replica symmetry breaking (RSB) schemes [25, 26] is a challenging and interesting 
future work. 

This work was partially supported by Grants-in-Aid from the MEXT, Japan, Nos. 
13680400, 13780208 and 14084206. 
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(a) (b) 




S[j Si 2 Si 3 Si 



FIG. 1: (a): Graphical expression of SG models in the case of K = 3 and C = 4. In this expression, 
each spin Si denoted as Q is linked to C = 4 couplings (□), each of which is connected to K = 3 
spins. and represent sets of indices of spins and couplings that are related to J M and Si, 

respectively. In the figure, £(/u) = {^1,^2^3} an d = {fii, /J,2, /X3, /X4}. (b): Cycles in a graph. 

A cycle is composed of multiple paths to link an identical pair of nodes. It is shown that BP can 
provide the exact spin averages in a practical time scale if a given graph is free from cycles || . 
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(a) 



(b) 



FIG. 2: Time evolution of macroscopic variables M* = (1/iV) £/=i mf and Q* = (1/iV) Ezli(™f) 2 
in the SK model for the BP updates (0) and (Q) (BP: 0)> the reduced dynamics @ and (|To|) 
(RS: lines) and the naive iteration of the TAP equation ( pd[ ) (TAP: +) for (a) Jo = 1.5 and 
(b) Jo = 0.5 keeping J = 1 and T = 0.5. TAP is plotted only for Q t in the case of Jo = 0.5 
in order to save space. Each marker is obtained from 100 experiments for N = 1000 systems. 
The AT stability is satisfied for Jo = 1.5 but broken for Jo = 0.5. Irrespectively of the AT 
stability, the behavior of the macroscopic variables in the BP dynamics can be well captured by the 
reduced dynamics while the naive iteration of the TAP equation does not exhibit any convergence 
even in the macroscopic scale. Insets: Squared deviation of spin averages between the successive 
updates D t = (l/N)Y$Li(mj — m\ x ) 2 is plotted for the BP dynamics. The deviation vanishes 
to zero indicating convergence to a fixed point solution for Jo = 1.5 while remains finite signalling 
instability of the fixed point for Jq = 0.5. 
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FIG. 3: Phase diagram for K = 2 and C = 4 suggested by eqs. (15) and (|i~6|). P, F and 
SG stand for the paramagnetic, the ferromagnetic and the spin glass phases, respectively. The 
boundary between F and SG is just a conjecture. In order to examine the validity of this diagram, 
100 Monte Carlo experiments were performed for N = 2000 systems at conditions denoted by +. 
For each condition, frequencies of macroscopic magnetizations M a = (1 / TV) Yli=i m f an d overlaps 
Qab = (l/-^0 Ya=i r nf' m \ ( a > ^) were evaluated, where mf is the average of Si obtained from 
20000 Monte Carlo steps per spin for experiments a,b = 1,2,..., 100 (insets). For both of the two 
conditions in P, all of M a and Q a b fall into the first bin. On the other hand, sharp peaks indicate 
the order to the ferromagnetic state in F and a broad distribution of Q a b signals the breaking of 
the replica symmetry in SG. 
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